Diffusion and binding of finite-size particles in confined geometries 
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Describing the difFusion of particles through crowded, confined environments with which they 
can interact is of considerable biological and technological interest. Under conditions where the 
confinement dimensions become comparable to the particle dimensions, steric interactions between 
particles, as well as particle-wall interactions, will play a crucial role in determining transport 
properties. To elucidate the effects of these interactions on particle transport, we consider the 
diffusion and binding of finite-size particles within a channel whose diameter is comparable to the 
size of the particles. Using a simple lattice model of this process, we calculate the steady-state 
current and density profiles of both bound and free particles in the channel. We show that the 
system can exhibit qualitatively different behavior depending on the ratio of the channel width to 
the particle size. We also perform simulations of this system, and find excellent agreement with our 
analytic results. 

PACS numbers: 02.50.-r, 05.60.-k, 05.40.-a, 87.15.hj 



I. INTRODUCTION 

Recent advances in technology and a burgeoning inter- 
est in biological systems have generated a great deal of 
interest in understanding particle transport in crowded 
environments. There are many biological processes where 
the diffusion of particles through the crowded environ- 
ment of the cell is important. Examples of these in- 
clude: The transport of material through ion chan- 
nels P, 0], mitochondrial and bacterial porins and 
nuclear pores [3| ; and the diffusion of enzymes and macro- 
molecules through microtubules and microtubule bun- 
dles [E, Q . In addition, cells themselves diffuse through 
confined environments, such as the movement of red 
blood cells and leukocytes through small blood vessels 0] ■ 
Finally, the transport of organelles between cells has been 
shown to occur through narrow "nanotubes" that con- 
nect the cells M. In technology, microfluidic devices and 
techniques (9l.Tl0j have immense potential for a variety 
of applications, including miniature biological assays for 
diagnostics and basic research Other applica- 

tions where such considerations would be important in- 
clude diffusion through carbon nanotube s Il3l and mi- 
croporous materials such as zeolites 0, [3 Il5|. as well as 
the diffusion of colloidal particles through narrow chan- 
nels 

To date, much of the theoretical effort on particle diffu- 
sion in confined environments has focused on two extreme 
limits. The first limit is the "single-file diffusion" case, 
where the effects of confinement are so severe that steric 
interactions between the particles prevent them from dif- 
fusing past one another [3, [H, [l^ HO) HH ■ In many sys- 
tems, however, steric interactions are important, but the 
level of confinement is not so extreme, allowing the diffu- 
sion of a small number of particles past one another. At- 
tempts have been made to model such systems using both 



continuum [24I and lattice-based approaches [23|, |2J, |25| , 
but these studies address regimes close to the single-file 
limit using either perturbative, quasi-single-file models 
or two-file lattice models. In the opposite limit, steric 
interactions between the particles are ignored, but other 
effects of the confining environment are taken into ac- 
count. For example, the friction between the particles 
and the confining walls can lead to hinderance of diffu- 
sion Q. Also, the variation of the cross-sectional area of 
a channel has been shown to lead to a generalized one- 
dimensional diffusion equation known as the Fick- Jacobs 
equation [2^ [23, [H, [1^, [s^]- Therefore, it is clear that 
understanding the general problem of confined diffusion, 
where the degree of confinement is lessened but steric in- 
teractions remain important, remains an important and 
under-explored endeavor. 

Another important effect in these confined systems 
that has not received extensive theoretical attention is 
the effect of specific and non-specific interactions be- 
tween the particles and the confining environment. For 
example, enzymes can bind to and chemically modify 
specific sites in tubulin when diffusing through micro- 
tubules @, i. Also, microfiuidic channel walls can be 
functionalized to allow for the binding of ligands in order 
to enable detection [HI, [l^l- Furthermore, electrostatic 
and van der Waals forcese are non-specific particle-wall 
interactions that, when combined with the effects of sur- 
face roughness and charge inhomogeneity of the walls, 
can lead to localized, transient binding of the particles 
to the walls. 

In this paper, we develop a simple model to study the 
diffusion of finite size particles through narrow channels 
with functionalized walls to which the particles can re- 
vcrsibly bind. We consider the limit in which the channel 
width is a few times larger than the particle dimensions, 
so that the particles can diffuse past each other relatively 



easily. On the other hand, the binding of particles on the 
channel walls can cause a bottleneck, effectively narrow- 
ing the dimensions of the channel for unbound particles. 
In section ini we describe in detail the setup of our basic 
model and (in conjunction with Appendix A) the ana- 
lytical procedure utilized to solve it. We also describe 
the simulations used to test our theoretical predictions. 
In section IIIII we describe the simplest case of diffusion 
in the absence of reversible binding, and make connec- 
tions to both the standard results for bulk diffusion and 
to the diffusion of particle through a channel of varying 
cross section [2^, [23, HO] ■ We then consider the 

effects of reversible binding on particle transport through 
the channel. We discuss two different cases, depending 
on the diameter of the channel relative to the particle 
size. In section ITVl we consider the case where the chan- 
nel diameter is small enough that it can be completely 
blocked by bound particles. We analyze the flux of parti- 
cles through the channel and the densities of bound and 
unbound particles within the channel, and show that our 
simulation results are in good agreement with the analyt- 
ical predictions. We also show that corrections to mean- 
field theory are necessary to account for the observed 
transport properties. In section |V] we consider the case 
where the channel is too wide to be completely blocked 
by bound particles. We show that the transport proper- 
ties, which in this case are adequately described by the 
mean-field theory, are significantly modified relative to 
bulk diffusion and to the case considered in section [TVl 
We conclude with a discussion of our results and their 
implications for biological and technical applications. 



II. MICROSCOPIC MODEL 

Consider a channel of cross-sectional area at in which 
particles of diameter S can both diffuse and bind to the 
channel walls. We partition the cross section into a num- 
ber of bins labeled by the index j, as illustrated schemat- 
ically in Fig.[TJ The ratio 4at/(7r(5^) determines the max- 
imum number N of particles that can fit within any given 
cross section of the channel; in our model, N is the num- 
ber of rows in the channel. We label sites along the 
axis of the channel with an integer i = 1, ...,7Vl, where 
Nl = L/S, L being the length of the channel. The in- 
dex j ranges from I, ...,N. Due to the steric interactions 
between particles, each site can be occupied by at 
most one particle. If a row in the model corresponds to a 
region adjacent to the walls of the channel, the particles 
can reversibly bind to the sites in this row. By varying 
the ratio w/S, where w cx y/Ht is the width of the chan- 
nel, we can see that there are two distinct cases we need 
to consider. When w/S ~ 1, every row j is an "exte- 
rior" row that lies along the wall of the channel. In this 
case, which we call the "no-hole case," it is impossible for 
particles to diffuse through cross sections of the channel 
which contain the maximum number of bound particles. 
When w/S 1, however, there are a certain number Nh 



j=N 




FIG. 1: Schematic illustration of the interior of the channel 
for a typical particle distribution of bound particles (shaded 
sites), unbound particles (cross-hatched sites), and unoccu- 
pied (white) sites. 



of "interior" rows where particles cannot bind. In this 
case, which we call the "hole case," particles can always 
freely diffuse in the center of the channel, even in regions 
where the maximum number of particles are bound to 
the channel walls. 

For both of the scenarios described above, the system 
evolves forward in time as a Markov process. That is, 
in a discrete time step At, each particle stochastically 
determines which (if any) of its possible moves it will 
attempt using a given state-independent probability for 
each move. Due to steric interactions, however, the cho- 
sen move can occur only if the particle attempts to move 
to an unoccupied site. If the site is occupied, the parti- 
cle does not move in that time step. There are several 
possible moves that each particle can attempt in a sin- 
gle time step. In general, the diffusion of particles in a 
channel can be "asymmetric" [2, 18, 19, 20, 21], with dif- 
ferent rates for hopping to the left and right. Also, the 
diffusion constant can inprinciple vary with the distance 
from the channel walls (i.e. it can depend on the row 
index j). In this paper, however, we assume that the dif- 
fusive landscape for the unbound particles in the channel 
is completely flat. In other words, we consider the sym- 
metric diffusion process exclusively. Furthermore, this 
assumption implies a symmetry between the rows in the 
channel that requires the left and right hopping rates to 
be independent of the row index j. Thus, an unbound 
particle at any site (i, j) can hop to the site (i ± l,j) 
with probability phop- In addition, an unbound parti- 
cle at (i,j) can hop to the site (i,j') with probability 
Phopiji j')- In principle, the functional dependence of 
Phop(j, i') on the rows j, j' must encapsulate the geome- 
try of the channel. The determination of these rates for 
arbitrary channel geometries could prove difficult, though 
the aforementioned symmetry between the rows requires 
Phop(j, j') = Phopij',.])- Fortunately, we shall see that a 
detailed knowledge of these rates will not be necessary 
to find the quantities of interest in this paper. Finally, if 
the index j labels an exterior row, an unbound particle 
at can bind to that site with probability pon, while 
a bound particle can unbind with probability poff ■ 

To determine the particle profiles, we need to spec- 
ify the boundary conditions at either end of the channel. 
Throughout this paper, we place the left end of the chan- 
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nel in contact with a bath of particles in the bulk, which 
sets the number of particles at z = 1. At the right end 
of the channel, we place an absorbing boundary (i.e. an 
infinitely dilute bath), which forces the particle profiles 
to vanish ai i — N^. 

Analytically, the time evolution of this system can be 
described by the master equation: 

i^p {t + At)) - \m) ^ \^m) = K m)) , (i) 

where K is the evolution operator and is the "state 

vector" of the system at time t: 

m))^ E E (2) 

{Mi_j}=0,l {bi_j} = 0,l 

Here, u^j is the number of unbound (bound) par- 

ticles at site (i,j), the "eigenstate" \s) = ^ij}) 
enumerates one specific configuration for all of the sites 
in the channel, and P[s]t] is the probability of finding 
the system in the eigenstate \s) at time t. The deriva- 
tion of the evolution operator K is given in Appendix [XI 
We note that although the sum in Eq. ^ includes physi- 
cally unallowable states (specifically, states with multiple 
particles occupying the same site), we can always force 
the probability of such states to vanish at all times (see 
Appendix [X| . 

Using the evolution operator K , we can compute the 
time evolution of the expectation value (i.e. the thermal 
average) of any physical observable. In this paper, we 
will focus excusively on the variation of physical observ- 
ables along the channel axis, rather than their variation 
within a given cross section of the channel. In particu- 
lar, we are interested in the total number of bound and 
unbound particles in the channel as a function of the (di- 
mensionless) distance i along the channel axis, as well 
as the current, i.e. the total number of particles per unit 
time passing through a given cross section of the channel. 
These quantities can be related to the expectation values 
of two operators, Cb{i,j) and Cu{i,j)- Cb{i,j) gives 1 if 
the site (i, j) is occupied by a bound particle, and oth- 
erwise, while Cu{i,j) gives 1 if the site {i,j) is occupied 
by a unbound particle, and otherwise. The evolution 
equations for these operators are derived in Appendix [XI 

In order to test the validity of our analytic results, 
we have performed simulations of the lattice model de- 
scribed above. All simulations were done on a rectangu- 
lar lattice like the one illustrated in Fig. [l] with N — 5 
and = 100. For simplicity, we set all of the lat- 
eral and longitudinal hopping probabilities to be equal, 
PhopU^j') = Phop- For every time step in the simulation, 
each particle in the channel is visited once, in a random 
order. During each particle visit, it is first determined 
what (if any) move that particle will attempt, using the 
rules and probabilities defined above. If that move is 
allowed - that is, if it does not lead to any multiply oc- 
cupied sites in the lattice - then it is performed; if that 
move is not allowed, then the attempt fails. This process 
is then repeated for the next (randomly chosen) particle. 



until all of the particles have been visited once during the 
time step. 

The boundary conditions in the simulation are set as 
follows: First, particles that leave either end of the chan- 
nel do not return. This alone sets the absorbing bound- 
ary condition at the right end of the channel. To set 
the boundary condition at the left end of the channel, 
we need an influx of particles into the channel at that 
end. To provide this influx, we stochastically attempt - 
with probability Pontor - to insert a single additional par- 
ticle into the left end of the channel at the end of each 
time step. If it is determined that an attempt should be 
made, then one of the N sites in the column i = 1 is cho- 
sen at random as the particle entry point. If that site is 
empty then the new particle is added there; if the site is 
occupied then the attempt fails. The value of the prob- 
ability Pcntcr sets the number of particles in the column 
i — 1, which must be measured in order to compare the 
simulation results to the theoretical solutions. 



III. DIFFUSION WITHOUT BINDING 

Before considering the full problem of diffusion and 
binding of finite-size particles inside a channel, we first 
consider the limit PomPoS ^ 0, in which the diffusing 
particles cannot reversibly bind to the channel walls. It 
is possible, however, to have an initial, stationary distri- 
bution of irreversibly bound particles in this limit. 

We first consider an initial condition with no bound 
particles. In this case, the problem reduces to the dif- 
fusion of finite-size particles with excluded volume inter- 
actions through a channel with a uniform cross section. 
Here, the only quantity of interest is the number of par- 
ticles along the channel, Nu{i,t) = {Cu{i,j))- Using 
the results of Appendix [XI it is straightforward to show 
that, in the continuum limit (At, S ^ 0) Eq. (|A10[) yields 
the standard diffusion equation for phantom (i.e. point- 
like) particles with: 

dtXuix,t) = DX'l{x,t), (3) 

where x = iS is the continuous position along the channel 
and D — lira^t,s^o PhopS^ / At is the diffusion constant. 
Here, Xuix^t) = liuis^o Nu{i,t) /S is the number of par- 
ticles per unit length at position x; that is, A„(a;, t) dx is 
the number of particles between x and x + dx. 

The fact that excluded volume interactions do not al- 
ter the simple diffusive behavior Eq. ^ of the particle 
profile is due to the assumed symmetry of the particle dif- 
fusion constant along the channel axis. Indeed, it is well 
known that excluded volume interactions do not affect 
the bulk diffusion equation when the diffusion constant 
is independent of position, even in the single-file limit 
(i.e. the symmetric exclusion process) In the case 

of an asymmetric exclusion process, where the hopping 
rate from i to i -I- 1 is different from the hopping rate from 
i to i — 1, excluded volume interactions do indeed affect 
the bulk diffusion equation ^3^. Furthermore, excluded 
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volume interactions do play a role in the behavior of in- 
dividual particles in the channel (i.e. tracer diffusion), 
even for symmetric diffusion processes [3l|. Finally, we 
note that the terms in Eq. (|A10[) for hopping within a 
given column i - i.e. the terms oc Phop(j,i') ~ cancel ex- 
actly in Eq. This occurs for arbitrary values of the 
lateral hopping rates PhopO, j'): long as these rates 
are symmetric, PhopO, j') = PhopW , j)- Thus, the lateral 
diffusion of the particles within the channel has no effect 
on the particle profile along the channel axis. 

We can also use the limit PomPoS ^ to study the dif- 
fusion of particles through a channel with a cross section 
that varies on length scales much longer than the parti- 
cle size. To do so, we choose initial conditions such that 
{Ch{i,j)) = ni,{i,j) is a fixed function that represents 
the varying cross section of the channel. Then Eq. (jA10[) 
becomes 

ANuit,t)^Pi,,pY,Y,[{nuii±l,j)-nuit,m (4) 
± 3 

X [1 - nb{i±l,j)][l - UbiiJ)]. 

If we assume that the cross section varies on length 
scales much longer than the channel radius, we can ap- 
proximate the distribution of particles within a given 
cross section of the channel by a uniform distribution. 
This is known as the "local equilibrium approxima- 
tion" "28]. In this hmit. 



Nu{^) 



[N-m{i)]' 



(5) 



where Nt{i) = Y.jMhj)- Using Eq. ([5]), Eq. Q be- 
comes 



ANu{i,t) =PhopX! 



A{i±l) A{i) 



g(z±l), (6) 



where the "cross-sectional area" A{i) = N — Ni,{i) is 
the number of particles that can occupy row i, and the 
"permeability" Q{i±l) = Y.j['^-nb{i±l,jm-nb{i,j)]. 
Now, the assumption that the cross section of the channel 
varies on length scales much longer than the particle size 
d implies that the function nb{i,j) is a slowly varying 
function of the row index i. Therefore, in the continuum 
limit nb{i ± 1, j) = nb{i,j) + 0{5), and the permeability 
is given by 

Q{i±l)^Y.^l-nb{i,3)f +0{5) 
j 

= ^ [1 - nb{i,j)] + 0{5) = A{i) + 0{5). (7) 

The second equality follows from the fact that 1 — 
= 0,1 always. Then in the continuum limit 
Eq. ^ becomes 



d\u{x,t) 
di 



-J'[x,t], 



(8) 



where the prime indicates a derivative with respect to 
the spatial variable x, and the current 



J(x,t) = -D { A{x) 



d_ 

dx 



^u{x,t) 

A{x) 



(9) 



Here, ^(a:) = lim^^o ^(*)/<5 is the continuum cross- 
sectional area; that is, A{x)dx is the maximum number 
of particles that can simultaneously occupy the region in 
the channel between x and x + dx. We note that the 0{5) 
terms of the permeability Eq. ([7]) vanish in the continuum 
limit 5 —t Q. 

Eqs. ([S]) and ^ are known as the Fick-Jacobs equa- 
tion, and have already been derived from the continuum 
diffusion equation for a cylindrically symmetric chan- 
nel [H, [13, m, m, iOl. Our result generalizes the va- 
lidity of the Fick-Jacobs equation to any channel with 
cross-sectional area A{x). In particular, a channel with 
changing area due to a change in the shape of the cross 
section will also exhibit Fick-Jacobs behavior if the shape 
change occurs slowly enough. 



IV. NO-HOLE CASE 

We now turn to the case in which particles can re- 
versibly bind to the walls of the channel. In this sec- 
tion, we consider the "no-hole" case, in which the diffu- 
sion of unbound particles through a region of the chan- 
nel can be completely blocked by bound particles in 
that region. In the language of the lattice model, par- 
ticles can bind to every site {i,j) of the channel. The 
state of this system can be described by two functions, 
Naii) = J2j {^a{i,j)), which give the expected number 
of bound {a — b) and unbound {a = u) particles at po- 
sition i along the channel axis. If we take the continuum 
limit, the evolution equations (|A6p and (jA10|) become, 
respectively, 

d 

— Xb{x,t) = konK{x,t) - kosXb{x,t), (10) 
ot 



d_ 

dt 



Xu{x,t) = ~konK{x,t) + kosXb{x,t) - J'{x,t), (11) 



where fcon.off = Pon,off/Ai and Xa{x,t) = 
lims^Q Na{i,t)/S. The discrete longitudinal current 
J(i, t) is given by 



J(i,t) 



Phop 

At 



Nuit + 1) - Nuii) 



(12) 



+ ( {c'uii,j)Cbii + - (c„(z + 1, j)a(z, j)> ) 



The continuous current in Eq. (|lip is related to the dis- 
crete current by J{x, t) — lim^^o Jih t)- Like the simple 
diffusion case, the terms for hopping within a given col- 
umn i cancel for arbitrary values of the hopping rates 
Phop{j,j') as long as Phop{j,f) = PhopifJ)- The first 
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two terms of Eq. (fT2l) give the lattice version of the cur- 
rent for phantom particles (i.e. Pick's Law), while the 
final two terms give the correction due to the fact that 
an unbound particle at j) cannot hop to a site (i± 1, j) 
if it is occupied by a bound particle. 

As mentioned in Section [Tll our model assumes that 
the diffusive landscape for the unbound particles in the 
channel is completely flat. In the no-hole case, this im- 
plies that any two rows of the channel are interchange- 
able. As a result, the steady state expectation value of 
any operator that acts on a single row j will be inde- 
pendent of j, as long as the boundary conditions arc in- 
dependent of j. In more detail, consider the probability 
P[s; t] of a particular configuration s of bound, occupied, 
and unoccupied particles in the channel arising at time 
t. After interchanging two rows j and j', we arrive at a 
configuration s' of particles that can arise with probabil- 
ity P[s'; t] at time t. Symmetry between the rows implies 
that P[s;t] ~ P[s';t]. When this symmetry is present. 



and the steady state current is 



(13) 



That is, averages involving operators of a single row j 
must be independent of j, since those averages involve 
summing over the probability distributions P[s;t] de- 
scribed above. Thus, we obtain the same expressions 
for the expectation value of any operator {0{j)) after 
interchanging j with j', implying Eq. p3p directly. 

In order to solve Eq. (Ilip , we must postulate a form for 
the two-point correlation functions appearing in Eq. (jl2p . 
The simplest form for these correlation functions is given 
by the mean-field approximation, in which the correla- 
tions between the two operators are neglected: 

{C,Xh])CS',j)) = (C„(*)) {C,(i')) ■ (14) 

Then the mean- field current, in the continuum limit, be- 
comes 



-DX'Jx,t), 



(15) 



where A = lim^^o N/ S is the maximum number of parti- 
cles (both bound and unbound) that can fit in the chan- 
nel, per unit length, 

At steady state, dtXu{x,t) = dtXb{x,t) — and the 
solution to Eq. (fTU]) is 



Xu{x) = KdXb{x), 



(16) 



where the disassociation constant = fcoff/fcon- Note 
that this solution is independent of the mean-field ap- 
proximation, Eq. p^ . Using Eqs. (fTSl) . (fTB|) . and the 

boundary conditions A„(L) = Xb{L) — 0, Ab(0) = X'^\ 
the steady state solution to Eq. pT|) is a simple hnear 
profile: 



xr{x) = x, 



(0) 



(17) 



J„,f(x) = -DX'^ix) 



xi:>^K,D 

L 



(18) 



Figures [2] and [3] show the bound particle steady-state 
profile Xb{x) and the the steady-state current Jq, respec- 
tively. In both figures, we can see distinct deviations 
of the simulations (points) from the predicted mean- 
field predictions (dashed lines). These deviations arise 
from the mean-field treatment of correlation functions 
{Cu{i, i)Cb{i' T ])) ■ In order to understand the physical 
processes that cause these deviations, let us first consider 
a quenched distribution of bound particles where the un- 
bound particle density has reached a steady state. In this 
case, there will be an absence of correlations between the 
bound and unbound particle distributions, because the 
bound particle distribution is invariant in time. How- 
ever, if we now consider a single binding or unbinding 
event, it is clear that there will be a transient change in 
the surrounding unbound particle density as it relaxes 
toward a new steady state distribution. Since particles 
bind and unbind on finite timescales, every such event 
will lead to a transient deviation in the correlation func- 
tions from their mean field value. This deviation will 
be particularly significant if a binding (unbinding) event 
blocks (unblocks) the entire cross section of the channel. 

To quantify the deviations from mean-field theory, we 
need to construct a dimensionless quantity that involves 
the two-point correlation functions appearing in Eq. (jl2p . 
Although the operators Ca{i,j) are dimensionless in their 
discrete form, their expectation values in the contin- 
uum limit Xa{x) have units of {length) ^. Therefore, 
we characterize the deviations of the expectation value 
(Cti(i± 1, j)Cb{i, j)) from its mean-field value with the 
dimensionless quantity 



_ {Cujt ± l,J)Cb{^,J)) - {Cu{i ± 1,3)) {Cb{t,j)) 
{Cu{i±l,j)){Cb{i,j)) 

(19) 

Here, we have indicated the independence of the func- 
tion x± from the row index j, which, as argued above, 
is known by symmetry. Physically we anticipate that 
X± will be a function of the dimensionless density of the 
bound and unbound particles. Assuming the deviations 
are small, we may expand X±(*) ™ powers of Nb/N and 
Nu/N and find, to lowest order. 



X±{i) = ei 



Nbii) , Nuit±l) 



N 



£2- 



N 



(20) 



Here, ei and £2 are parameters that encapsulate the de- 
gree of deviation from mean-field behavior. Note that 
these parameters are not universal: in principle, they de- 
pend on the various hopping rates. 

Given Eq. ([^D|) . the current can be written as J{x, t) = 
Jtn{{x,t) + Jcori{x,t), whcrc Jni{{x,t) is givcu by Eq. ([H 
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FIG. 2: Dimensionless bound particle density Ai,(a:)/A for 
the no- hole case obtained from the simulations (points), as 
compared to the theoretical predictions both at the mean- 
field level (dashed lines) and including the correlations (solid 
lines). For all data shown, N — 5, Pon = 0.01, Pos = 0.001, 
andphop = 1/150; for the simulations, we also set Phop(j, j') ~ 
Phop for all The value of pcntcr for each simulation sets 

the left-end boundary condition A^""* for the theoretical curves. 
The values of Pontcr and A^^', respectively, are: 0.006, 0.60 
{top, squares); 0.02, 0.79 (top, triangles); 0.01, 0.70 (bottom, 
circles); and 0.5, 0.90 (bottom, diamonds). For the theoretical 
predictions that include the effects of correlations, all of the 
curves use the same fitting parameter, e — 0.27. 



Dei 
A2 



[Xuix,t)d,Xl(x,t) ~ Xl(x,t)X'^(x,t)] 



and 

t-^corr {x , 

D( 

T2 



The steady-state relation given by Eq. (jT6|) still holds, 
since (as noted above) it is independent of the mean-field 
approximation. However, the steady-state bound particle 
profile and current are now given by, respectively, 



+ ^ [>^b{x, t)d,Xl(x, t) - Xl(x, t)X',(x, t)] . (21) 



X,(x)^XT(x) 





2 











(22) 



0.004 




0.2 0.4 0.6 0.8 



FIG. 3: Dimensionless steady-state current Jo/fchop obtained 
from the simulations (points) for various values of the left-end 
boundary condition X^^^ /K, as compared to the theoretical 
predictions both at the mean-field level (dashed lines) and 
including the correlations (solid lines). All parameter values 
are identical to those used in Fig [5] 



and 



J(x,t) = Jo = 



DKaX 



(0) 



L 



1 



A 



(23) 



where e = (— ei + ii'de2)/3. Here, we have only retained 
the terms linear in e in Eq. ((22|) , in order to be consistent 
with the expansion of the correlation functions Eq. (|20p . 
Thus, we can see that our proposed form of the deviations 
from mean-field behavior has only one fitting parameter, 
e. 

Figures [2] and [3] show the resultant fits (solid lines) of 
Eqs. ([22]) and ((23)) . respectively, to the simulation data . 
In order to determine the value of the fitting parameter 
e, we fit Eq. ([23|) to the current data shown in Fig. [3l 
and then use this value for all of the particle profiles 
shown in Fig. [2l We can see that this gives excellent fits 
for all of the profiles and for the current measured by 
the simulations. Thus, our postulate for the form of the 
deviations from mean-field theory captures the deviations 
seen in the simulations using only a single parameter fit. 



V. HOLE CASE 

We now turn to the case in which the ratio of the 
channel diameter to the particle size is large enough 
that particles can always diffuse through a "hole" in 
the center of the channel, even when all of the bind- 
ing sites in a given cross section of the channel are oc- 
cupied by other bound particles. In our lattice model, 
this corresponds to two distinct types of rows: exterior 
rows, in which particles can diffuse and reversibly bind, 
and interior rows, in which the particles can only dif- 
fuse. Thus, the state of the system can be described by 
three operators: Na{i) = {Ca{h j)) 8j for a = u,b 
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gives the expected number of bound (a = b) or un- 
bound {a = u) particles in the exterior rows of column 
i; Nh{i) — J2j {C^uihj}) (1 — ^j) gives the expected num- 
ber of unbound particles in the interior rows of column i. 
Here, 9j — 1 for the exterior rows, and for the interior 
rows. 

For both the case of diffusion without reversible bind- 
ing and the no-hole case considered above, the lateral 
hopping terms in the evolution equations of interest can- 
celled one another exactly. This cancellation occurred 
for arbitrary values of the lateral hopping probabilities 
PhopOij')) subject only to the symmetry requirement 
Phop(j, j') = Phop{j',j)- As we shall see below, however, 
the distinction between the interior and exterior rows in 
the hole case causes some of the lateral hopping terms 
- specifically, those terms corresponding to the diffusion 
of particles from the interior rows to the exterior rows 
(and vice-versa) - to remain in the relevant evolution 
equations. Consequently, we need to make a simplifying 
assumption about these hopping probabilities. Since the 
primary focus of this paper is the effects of steric interac- 
tions on the diffusion of particles in confined geometries, 
we need not consider cases where the number TVh of in- 
terior rows is large; the particle diffusion in such systems 
can be adequately described by the standard diffusion 
equation for phantom particles. Geometrically speaking, 
when iVjj is small, all of the exterior rows will be approx- 
imately equidistant from any given interior row. There- 
fore, we will assume that the probability of hopping from 
any exterior row to any interior row (and vice-versa) is 
given by phop = const. This assumption does not apply 
to the hopping probabilities from one interior row to an- 
other interior row, or from one exterior row to another 



exterior row: these rates remain arbitrary, except for the 
usual the symmetry constraint PhopU^-f) — Phop(j', j)- 
Thus, we assume that the lateral hopping rates are of 
the form 



Phop{j,j') = Phop [1 - 6ij'] + (j>{jj')Qjj' 



(24) 



where (j){j,j') = 4>{j',j) and Qjj' — if 6j ^ 9ji - that is, 
if one row is an exterior row and the other is an interior 
row - and 1 if 9j = 9ji . 

If we take the continuum limit, the evolution equa- 
tion for the bound particle profile reduces to the same 
equation obtained in the no-hole case, Eq. (fTU]) . We can 
use the evolution equation for {Cu{i,j)), Eq. (jAlOp . to 
obtain the evolution equations for both of the profiles 
Nu{i,t) and Nii{i,t). Specifically, if we use the assump- 
tion Eq. for the lateral hopping rates, it is straight- 
forward to show the terms oc cancel one another 
exactly in both evolution equations, so that 



NS) 
N -Nh 



Nh{i) 



N -N, 



H 



H 



(1 - 0,) . 



(25) 
(26) 



In contrast to our discussion of the no-hole case, we can- 
not exchange an interior row with an exterior row, al- 
though we can interchange interior (and exterior) rows 
amongst themselves. This latter symmetry ensures that 
the expectation values take on only one value on all 
interior rows and another on all exterior rows. Using 
Eqs. ([24| and (|25|) . it is straightforward to show that 
Eq. (jAlOp gives the mean-field evolution equations 



d_ 

dt 



\u{x, t) = -konK{x, t) + /coff Ah(x, t) - DA [\u{x, <)Ah - Xh{x, t) {A-Ah- Xb{x, t))] - j;(x, t), (27) 



— A?, (a:, t) = DA [Xu{x, t)AH - \h{x, t) (A - - Afc(.T, t))] - J'^{x, t), 



(28) 



where \a{x^ t) = lima^o Na{i, t)/6 for a = u,h, b, Ah — lim^^o Nh/S, D = limAt,5^o khopS'^ / (N At) , and the currents 

1 



D 



KX^,t) + 



A- A 



H 



{\u{x,t)X',{x,t) - Xt{x,t)X'^{x,t)) 



Jh{x,t) = -DX',^{x,t). (29) 



r 



At steady state, we can see that Eq. (IT6|) still holds, 
and that the total longitudinal current Jtot(2;) = Juix) + 
Jh{x) is constant: 



a = u,b,h), implies that 



Xhix) 



Jo 
D 



[L- x) - KdXb{x). 



(31) 



Jtot(x) = -D [X'Jx) + X'^ix)] = Jo. 



(■30"^ To solve for the steady-state current Jq, we need an ad- 

ditional boundary condition relating the particle profiles 
in the exterior and interior rows. Outside of the channel, 
This, along with the boundary condition Xa{L) — (for there is no net flux of particles in the direction perpen- 
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0.012 



FIG. 4: Dimensionless bound particle density Ai,(x)/A for the 
hole case obtained from the simulations (points) , as compared 
to the mean-field theoretical prediction (solid lines). For all 
data shown, Nh ~ 1; the remaining parameter values are 
identical to those used in Fig. [2] The values of pcntor, and the 
resultant values of aJ,''^ are, respectively: 0.004, 0.41 (circles); 
0.008, 0.53 (squares); 0.02, 0.63 (diamonds); and 0.5, 0.72 
(triangles). 



dicular to the axis of the channel. Since the current must 
be continuous across the boundaries of the channel, the 
lateral current at the channel ends must vanish. Specif- 
ically, the rate of particles hopping from the interior to 
the exterior rows must equal the rate of particles hop- 
ping from the exterior to the interior rows at the channel 
ends; that is, the term in brackets in Eq. (pS)) must vanish 
there. This is trivially satisfied at the right end at steady 
state, since all of the particle profiles vanish there. At the 
left end, we must have 



KdXb{0)^H = Aft(0)(A -Ah- Xb{0) 



(32) 



Such an equation is not necessary in the no-hole case, 
since in that case symmetry dictates that the lateral cur- 
rent vanishes everywhere. Using Eq. ([5^ . the steady 
state current is given by 



Jn = 



L 



A-Afc(O) 



A- A 



H 



XbiO) 



(33) 



The remaining ODE for the steady state profiles can 
be obtained by combining Eqs. (^5]) . and (PT|) : 



K^Ux) = DA 



(^(i -x)- K^Xbix)) [a -Ah- Xb{x)) 



KdXb{x)AH 



(34) 



This nonlinear ODE must be solved numerically using 
the boundary conditions Ab(0) = X^^^ and Xb{L) — 0. 

Figures [4] and \5\ show the simulation results (points) 
for several steady-state bound particle profiles A(,(x) and 
the the steady-state current Jq, respectively, for the hole 
case. We can see that the mean-field predictions (solid 
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FIG. 5: Dimensionless steady-state current Jo/fchop obtained 
from the simulations (points) for various values of the left-end 
boundary condition X^^-'/A, as compared to the mean-field 
theoretical predictions Eq. (|33|l [solid lines]. All parameter 
values are identical to those used in Fig|4l 



lines) for the particle profiles and the current show ex- 
cellent agreement with the simulation results with no fit- 
ting parameters. Thus, in contrast to the no-hole case, 
the effects of the correlations that are ignored in the 
mean-field approximation are negligible in the hole case. 
Specifically, we can see from original evolution equation, 
Eq. (jAlOp . that the relevant two-point correlation func- 
tions are of the form (C„(z, j)Ch(i', j')). In the no-hole 
case, the correlation of these two operators can be signif- 
icant, because the hopping of unbound particles can be 
completely prevented by a large number of nearby bound 
particles. This is never true in the hole case, however, be- 
cause unbound particles can always diff'use through any 
region of the channel using the interior sites. Thus, the 
correlations in the hole case are always negligible, and 
mean-field theory provides an excellent description of this 
system. 



VI. DISCUSSION AND CONCLUSION 

Diffusive transport driven by concentration gradients 
can occur under varying degrees of confinement, depend- 
ing on the relative size of the channel and the diffusing 
entity. In one extreme, the particle size is much smaller 
than its surrounding environment, and the transport be- 
havior obeys the standard Pick's law. In the other ex- 
treme, the degree of confinement is severe, forcing the 
particles to diffuse in a single file. In this paper, we 
have presented the first steps toward a systematic un- 
derstanding of diffusion in systems with an intermediate 
degree of confinement, where steric interactions prevent 
a large number of particles from diffusing freely past one 
another. For symmetric diffusion through a channel with 
a uniform cross section, we recovered the standard bulk 
diffusion equation, which neglects steric effects. This is 
not surprising: Although steric confinement does affect 
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the single particle dynamics of diffusing particles, it does 
not alter the bulk diffusive behavior for a symmetric dif- 
fusion process, even in the single-file limit [3l|. We also 
examined diffusion in a channel with a slowly varying 
cross section, and were able to derive a g eneralized form 
of the Fick- Jacobs equation [H 113, [IS Is^l for the 
diffusion of particles in a channel with a varying cross- 
sectional area. 

In such confined geometries, ubiquitous specific or non- 
specific interactions between the diffusing objects and the 
channel walls fundamentally alter the transport dynam- 
ics. When both steric confinement and reversible binding 
to the channel walls are present, there are two qualita- 
tively different cases that arise as the ratio of the channel 
width to the particle size is varied. In the first (no-hole) 
case, the cross section of the channel is wide enough to 
accommodate several particles, but still narrow enough 
that the diffusion of particles through a particular region 
of the channel can be completely blocked by bound parti- 
cles in that region. Our simulation results for the bound 
particle profile indicate a monotonic decrease from the 
proximal (left) to the distal (right) end of the channel. 
The steady state current increases as the density of bound 
particles at the left end increases, but begins to show a 
saturation behavior at high values, stemming from the 
fact that the system spends more time in configurations 
in which the channel is completely blocked. Interestingly, 
our simulation data for both the particle density profile 
and the total current at steady state reveal significant 
deviations from our mean field predictions, especially at 
higher overall particle densities. This is due to the fact 
that binding events that completely block the channel, as 
well as unbinding events that relieve this blockage, lead 
to significant deviations in two point density correlation 
functions from their mean field values. By taking these 
deviations into account by means of a single dimension- 
less parameter, we can reproduce both the particle den- 
sity profiles and the steady state currents seen in all of 
our simulations. Being able to characterize the parti- 
cle profile and current across a wide range of concentra- 
tion gradients across the channel with a single parameter 
is bound to be extremely useful in predicting transport 
behavior in systems where a limited amount of data is 
available. 

For the second (hole) case, the channel is wide enough 
to allow diffusion of particles through its center, even 
in regions where there is a saturating coverage of bound 
particles on the wall. Our simulation data in this case 
is in excellent agreement with the predictions from mean 
field theory. Since the channel cannot be blocked under 
any circumstance, the two-point correlations between the 
bound and unbound particles are much weaker than in 
the no-hole case, making the deviations of the particle 
profiles and current from their mean-field values negligi- 
ble. Both the particle profiles and the current show qual- 
itative differences from the no-hole case. At high values 
of the left end particle density, the bound particle profile 
shows a plateau phase near the left end before dropping 



to zero at the right end. This indicates that as the con- 
centration gradient across the channel is increased, the 
particles bind and effectively coat larger and larger re- 
gions of the channel wall, starting at the left end. This 
is also reflected in the steady state current, which shows 
a remarkable biphasic behavior as the left end particle 
density is increased. At low densities, the current shows 
only modest increases as the density is raised; at high 
densities, on the other hand, the current rises sharply 
for increasing densities. This is due to the aforemen- 
tioned coating effect of the bound particles at high den- 
sities, which forms a non-sticky layer along the channel 
walls. The diffusive behavior then becomes akin to that 
of phantom diffusion in a non-sticky channel (albeit of 
a smaller cross section), resulting in strong increases in 
the current at high concentrations. This kind of strong 
biphasic behavior could have important implications in 
a variety of biological systems, where it could be used as 
a regulatory or sensory mechanism. In artificial systems, 
one could potentially tune the system properties to gen- 
erate a desired strongly non-linear dependence between 
the current and concentration gradient. Finally, it is im- 
portant to note that one could go from the no-hole to 
the hole case with a small change in the channel diame- 
ter (on the order of a particle size). That such distinct 
transport regimes are separated by such small changes in 
the geometry of the channel could also have wide-ranging 
implications. 

Promising avenues for further research include extend- 
ing our approach to include a bias in the diffusion (i.e. 
asymmetric diffusion), which could naturally occur as a 
result of electric fields, hydrodynamic flows or even bi- 
ased motion of molecular motors. Interactions between 
the particles themselves could also yield significant new 
regimes. We hope that our work on these under explored 
systems, where the intermediate degree of confinement 
and the particle-wall interactions lead to novel and qual- 
itatively different transport behaviors, inspires further 
theoretical, computational and experimental research on 
these very rich systems. 
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APPENDIX A: THE MASTER EQUATION 

The evolution operator K defined in Eq. ^ can be 
expressed in terms of the creation and destruction oper- 
ators a^{i,j) and a^{i,j) for bound and unbound parti- 
cles, respectively: 

4(*,j)|0;6),,-|i;6),, 

<(*,j)|i;^>.,, -|0;fo)„- (Ai) 

«b"(*,j)k;i>.,, = k;0)„. (A2) 
atihJ) \u; l),j = W;0)^ j = 0, 

where |u; &),- ^ specifies a state of the site {i, j) with u un- 
bound and b bound particles. It is important to note that 
these operators can create the unphysical state |l,l)j^, 
in which the site is occupied by both a bound and 
unbound particle. Therefore, we must take care to con- 
struct the evolution operator K so that only transitions 
between physically allowable states can occur. This will 
ensure that, as long as the initial condition of the sys- 
tem is a physically allowable state, the system will never 
evolve into unphysical states. 

We can construct number operators from the creation 
and destruction operators, na(i,j) — a^{i, i)a'^(i, j), 
where a = u,b. In terms of these number opera- 
tors, the operators Cb{i,j) = ni,{i,j)[l — n„(i, j)] and 
Cti(t,j) = n-uii,j)[l — ni,{i,j)]. We also define the oper- 



ator Ce{i,j) = [1 - nb{i,j)][l - nu{i,3)\, which gives 1 if 
the site («, j) is empty and otherwise. 

Consider a single event, transforming the system from 
a state \soid] to a state |snetu)j that is allowed to occur in 
a time step At [e.g. the hopping of an unbound particle 
from (i, j) to (i -|- 1, j)]. Using Eqs. H]) and (O and the 
orthogonality of the state vectors, (s| s') ~ (5s,s', 

A^[s;i]= E P{^-A{s\K\l) (A3) 

{ft;'^}=0,l {f4-^}=o,i 

where for any function /(t), A/(i) = ]{t + M) - ]{t). 
For every possible transition jsoid) — > \snew), there must 
be two terms in K . Both of these terms should be pro- 
portional to P\soi4_\i\^ since the frequency of the transi- 
tion \soid) — > \snew) clcarly depends on the probability 
of finding the system in the initial state \soid)- The first 
term accounts for the increase in P[s„etu;^] due to this 
transition. This term should give a positive contribution 
to the RHS of Eq. (jA3|) for s = s„eu;- Then it is clear 
from Eq. (jASp that this term should be positive and con- 
tain creation and destruction operators that transform 
\soid) \snew)- The sccoud term in K for this transition 
accounts for the decrease in P\soid\ i] due to the transi- 
tion \soid) — > \snew)- This term should give a negative 
contribution to the RHS of Eq. (jA3p for s — Soid- We 
can see from Eq. (|A3|) that this term should be negative 
and contain only number operators, so that the non-zero 
term in the sum on the RHS of Eq. (|A3p is proportional 
to P[soid',t]. Using these rules, it is straightforward to 
write down the evolution operator K for the system de- 
scribed in section ini Writing K = J2i j ^jj ; 



Ki^j ^\^Pon[a^{i,j)au {hi) ~ j) (1 - j))] + Pos[% {i,j)a+{i,j) - (1 - nu{i, j)) ni,{i, j)]j0j 

+ E^hop(j,i') (1 - nb{i,j)) (1 - nb{i,j')) [a:^{i,j')a+{i,j) - nu{i,j') (1 - nu{ij))] (A4) 

± 

I 



where 9j = 1 if binding can occur in the jth row, and if 
binding cannot occur. As required, this evolution opera- 
tor satisfies the constraint that only physically allowable 
states evolve in time. 

We can now use the master equation Eq. ^ to com- 
pute the evolution equations for any given operator. 
From Eq. 12]), it is clear that the normalization of the 
state vector \ip{t)) is given by (U^ | ipit)) = !> where 
= J2s I*)- Therefore, the expectation value of any 
operator O^J is (0,j) = Oij lipit)). Using Eq. (P), the 
evolution equation for (Oij) is given by 

A{0.,)^{r\0,,Km))- (A5) 



To compute the RHS of this equation, we note that if 
Ki'j' contains no operators that act on the site (i, j), then 
(U^l OijKiiji \ip{t)) = 0. Using this fact, it is straightfor- 
ward to derive the desired evolution equations for the 
operators of interest. In particular, the evolution equa- 
tions for Cb{i,j) and Ce{i,j) are, respectively, 

A(C,(z,j)) = [Pon {Cu{i,j))-PoS {Cb{i,j))]e,, (A6) 

A{Ce{i,j)} + A{Cu{i,j)) + A{Cb{i,j))=0. (A7) 

Since every site must either contain a single particle 
(bound or unbound), or be unoccupied, and {Cb{i,j)) = 
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if 6j ~ 0, the solution to Eq. \K1\ is simple: 

j)> = 1 - j)) - j)) (A8) 

Indeed, for any operator Oi' y , 

{0,yC,{i,j)) = {O,,,, [1 - C„(z, j) - Ck{i,])6,]) . 

(A9) 



Using this result, it is straightforward to show that the 
final desired evolution equation, for the operator Cu{i^ j), 
is given by 



+ {Cu{t,j)Cb{t,f)) 0y] +Phop^ [{Cu{t ± 1, j)) - j)) - {Cu{t ± l,j)Cb{i,j)) 0j + {Cu{i.])CS ± 1, j)> 

(AlO) 

I 
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